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Abstract 

We present a new numerical technique to solve large-scale eigenvalue problems. It is based on the projection technique, 
used in strongly correlated quantum many-body systems, where first an effective approximate model of smaller complexity 
is constructed by projecting out high energy degrees of freedom and in turn solving the resulting model by some standard 
eigenvalue solver. 

CNJ Here we introduce a generalization of this idea, where both steps are performed numerically and which in contrast to 
^jfjthe standard projection technique converges in principle to the exact eigenvalues. This approach is not just applicable 
^ to eigenvalue problems encountered in many-body systems but also in other areas of research that result in large scale 

eigenvalue problems for matrices which have, roughly speaking, mostly a pronounced dominant diagonal part. We will 

present detailed studies of the approach guided by two many-body models. 
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1. Introduction 

The solution of large, sparse eigenvalue problems is an 
important task in engineering, mathematics, and physics, 
particularly in the field of strongly correlated quantum 
many-body systems, such as the high-Tc superconduc- 
tors, and more recently cold atoms on optical lattices and 
coupled light-matter systems. The treatment of strong 
correlations between particles leads to exponentially large 
eigenvalue problems as a function of the system size. 

To begin with, algebraic eigensolvers like the Lanczos 
and Arnoldi methods (TJ |2] play an important role in 
many body-physics as well as in many other research ar- 
eas, mainly due to the fact that they are widely applicable, 
simple, effective and they yield numerically exact results. 
Implementations of the methods can be obtained from the 
Internet, e.g. the well-known ARPACK[^ routine which in- 
cludes an implicitly restarting Arnoldi algorithm. On the 
downside, only comparably small systems can be treated 
by these techniques. For most problems in quantum many- 
body physics, the corresponding geometric sizes are much 
too small. 

A couple of sophisticated numerical methods have been 
developed for such systems in recent years. The density 
matrix renormalization group (DMRG) [3] is a powerful 
algorithm for the determination of ground state properties 
of large systems, which are not algebraically feasible. Re- 
cently the approach has been embedded into a wider math- 
ematical framework, the matrix product states (MPS, [4 ). 
The approach is limited to ID or quasi- ID systems. An- 
other method is the cluster perturbation theory (CPT, [5]), 
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which constructs approximations to the Green's function 
of infinite systems in 1, 2, and 3 dimensions by exact treat- 
ment of small clusters and combining them by perturba- 
tion theory. An extension of CPT represents the varia- 
tional cluster approach (VGA, [6 ) which improves the re- 
sults by introducing variational parameters. An approach 
without systematic errors is given by the family of quan- 
tum Monte-Garlo techniques (QMG, [7 ). QMG simula- 
tions, which are based on high dimensional random sam- 
ples of some suitable probability-density function have an 
statistical error which declines with the sample size. They 
are applicable to fairly large systems and to finite temper- 
atures. The drawback of these methods is the so-called 
sign-problem^ that shows up especially in fermionic mod- 
els and which makes certain models or parameter regimes 
inaccessible [8 . The methods discussed so far are par- 
ticularly tailored for quantum many-body problems and 
cannot be applied easily, if at all, to matrices of other ap- 
plications. 

The methods for strongly correlated quantum many- 
body systems rely on the assumption that the model has 
only a few rather local degrees of freedom. Real ab-initio 
models for strongly correlated quantum many-body sys- 
tems are out of reach in any case, and it is inevitable to 
reduce the complexity of the system upon describing the 
key physical properties by a few effective degrees of free- 
dom, like in the multi band Hubbard model [9]. In many 
cases, these models are still too complicated and cannot 
be solved reliably neither analytically nor numerically. For 
these cases it has been proven very useful to construct 
effective Hamiltonians like the Heisenberg- or t J-model. 
They are obtained via the projection technique 
upon integrating out such basis vectors which correspond 
to high energy excitations, or rather which have very large 
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diagonal matrix elements in a suitable basis. Of course, 
it is desirable that the quantitative results of the effective 
model are close to those of the original model. But more 
than that, it is the qualitative generic physics of strongly 
correlated fermions at low energies, which one wants to un- 
derstand. Since the original model itself is tailored for that 
purpose and is already a crude approximation of the under- 
lying ab-initio Hamiltonian, it suffices to have an effective 
model that still includes the key physical ingredients in 
order to describe competing effects of strongly correlated 
quantum many-body systems. The standard projection 
technique defines the space of dynamical variables (basis 
vectors) which are of crucial importance for the low lying 
eigenvalues, which are marked by small interaction energy 
(small diagonal matrix elements). The residual dynamical 
variables (basis vectors) are treated in second order per- 
turbation theory. If the model contains too complicated 
terms, such as the density assisted next-nearest neighbor 
hopping in the t J-model, they are omitted. The resulting 
model has a significantly reduced configuration space and 
is solved by one of the above mentioned techniques. 

Here we present a numerical scheme, that represents a 
threefold generalization of the projection technique: a) it 
combines the two steps of the construction of the effective 
model and exact diagonalization, b) no approximations are 
made as far as the effective model is concerned and c) it 
allows a systemic inclusion of higher order terms up to 
the convergence to the exact result. The projection step is 
based on the Schur complement and results in a non-linear 
eigenvalue problem, which is solved exactly. 

The numerical projection technique (NPT), that shall 
be discussed in this paper, has several advantages. First 
of all, it is not necessary to formulate an explicit effective 
Hamiltonian, which is not trivial in more complex models 
involving several bands and other degrees of freedom than 
just charge carriers [11 . And last but not least, NPT al- 
lows to systematically go to higher orders, which becomes 
necessary if the coupling strength is merely moderate. 

NPT is applicable to other large scale eigenvalue prob- 
lems as well. The underlying matrix has to be sparse and 
the diagonal elements, after suitable spectral shift, have 
to fulfill the following criterion: a few diagonal elements 
are zero or small and the vast majority of the diagonal 
elements is (much) greater than the sum of the respective 
off-diagonal elements. 

The method is demonstrated by application to two 
representative problems of strongly-correlated quantum 
many-body physics, the spinless fermion model with near- 
est neighbor repulsion and the drosophila of solid state 
physicists: the Hubbard model with a local Coulomb re- 
pulsion. These models are used for benchmark purposes 
only and it is not the goal of the present paper to provide 
a detailed discussion of the fascinating physics described 
by these models. 

The outline of the paper is as follows. In the second 
section we present the two benchmark models. The nu- 
merical projection technique is introduced in detail in the 



third section and analyzed in the ensuing sections. 

2. Model systems and basic idea 

The approach that we present below is generally appli- 
cable if the matrix, for which the lowest eigenvalues shall 
be determined, can be split into parts of increasing diago- 
nal dominance. What this means in detail will be clarified 
using two examples of the realm of many-body physics, 
the spinless fermion model with strong nearest neighbor 
interaction and the Hubbard model for fermions. Both 
models are tailored to study effects of strong correlations 
of electronic systems, such as the Mott-insulators [12 , high 
temperature superconductors [13 , manganites [14 , just to 
name a few of the very many novel materials with fasci- 
nating many-body effects. The Hamiltonian of the spinless 
fermion model reads 

H = —t ^ alcLj + V ^ fiifij , 

where a] (a • ) denotes the creation (annihilation) operator 
for fermions at site i. These operators have the common 
fermionic anti-commutator relations (for details see [H]). 
The operator hi = old- is the particle number operator for 
site i. The bracket (ij) indicates that the sum is restricted 
to nearest-neighbor sites Xi and Xj. There is a total of L 
lattice sites x^, which are placed on a simple cubic lattice 
in one dimension in the present work. 

In the case of the Hubbard model, the spin degree of 
freedom is also taken into account and the Hamiltonian 
reads 

H = -tY^ dl^dj^ ^UY^ hi^hj^ . 

{ij)a i 

The creation- (annihilation-) operators i^ia) obtain 
a spin index a with two possible orientations i- The 
operator hia- = ^L^icr particle number operator for 

spin a at site i and the operator for the particle density at 
site i is given by hi = hi^. 

The physical behavior of both models is determined by 
the relative strength of the hopping parameter t and the 
two interaction parameters /7, which stands for the on- 
site Coulomb repulsion between fermions of opposite spin, 
and V, which represents the repulsive nearest-neighbor in- 
teraction. In order to simplify the discussion we denote 
both parameters by V. In both cases the most interesting 
case is that of (almost) half-filling. That is, the number 
of particles in the system is half the maximum capacity, 
which is L in the spinless fermion model and 2L in the 
Hubbard model. Both models are particularly interesting 
in the strong coupling case, i.e. for V > \t\. 

For both models we use the occupation number basis in 
real space in which the interaction part is diagonal and has 
a value N^V^ where is either the number of occupied 
nearest-neighbor sites (Nnn) in case of the spinless fermion 
model, or the number of double-occupancies (N^) in case 
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Figure 1: Illustration of the Hamilton matrix of a V-ordered occupa- 
tion number basis. The blue parts contain the hopping terms inside 
a partition and the potential on the diagonal and the green parts 
contain the hopping terms between the partitions. 



of the Hubbard model. In strong coupling, states with 
increasing values of N,^ are decreasingly important. In 
the projection technique [9l [10] effective strong coupling 
models are derived, in which the configuration space is 
restricted to the sector = 0, i.e. Nnn = {N^ = 0) and 
the infiuence of the higher sectors are taken into account 
up to second order in t/V. A prominent example is the 
the spin- 1/2 antiferromagnetic Heisenberg model which is 
obtained in the half-filled case of the Hubbard model or 
the tJ-model, which is the generalization away from half- 
filling. Already in the t J-model terms are neglected which 
belong to the same order in t/V and to the same Nd-sectoi. 
In models with more bands and degrees of freedom, the 
derivation of an effective strong coupling model can be 
very demanding, like e.g. in the spin-orbital model for the 
manganites [11 . 

In this paper, we exploit this idea numerically and gen- 
eralize it in such a way, that the infiuence of high sectors 
is taken into account recursively without the need of leav- 
ing terms out that in standard projection technique would 
complicate the resulting effective Hamiltonian. For obvi- 
ous reasons we will refer to this approach as numerical 
projection technique (NPT). 



3. Numerical Projection Technique 

In order to exploit the NPT, we reorder the basis vectors 
according to A^,^. The corresponding Hamiltonian matrix 
has a natural block structure (see Fig. [T]) corresponding 
to the sectors with N,^ = 0, 1, 2, . . .. 

Note, that the hopping of a particle can only change the 
number by one. Therefore, only neighboring sectors are 
coupled and the matrix has a tridiagonal block-structure. 



3.1. Projection Step 

We start out with a general projection based idea for 
solving an eigenvalue problem of a 2 x 2 block matrix 



Aox= (^^^ 



A 



(1) 



The blocks correspond to some suitable partitions of the 
vector space under consideration. The sizes of the two 
partitions shall be denoted by pi and respectively. It 
is not necessary that the two partitions cover the entire 
vector space of the original problem. 

We assume that A and E are 'small' compared to 
such that the lowest eigenvalue of Aq is predominated by 
the corresponding eigenvalue of A and the perturbation 
due to E and B can be included perturbatively. In order to 
quantify this idea, we will map the infiuence of the second 
partition into the first partition by a Schur transformation, 
similar to the projection technique [H [To]. I.e. we multiply 
Eq. [ijfrom the left with the matrix 



-EB- 
I 



resulting in 



S 

E^ 



= A 



-EB- 
I 



with S = A — EB ^E^ being the Schur-complement. The 
second line of the eigenvalue equation yields 

X2 = -{B -\)-^E^xi . (2) 

Inserting X2 into the first line leads to 

Sxxi := {A - E{B - \)-^E^)xi = Xxi . (3) 



Note that this equation is in principle exact, no approx- 
imations have been made so far. The original problem 
(Eq. [T]) has thus been projected into the subspace of the 
first partition which is of smaller size, but at the prize 
of a non-linear eigenvalue problem. By the above proce- 
dure we obtain only those eigenvectors x = (xi,X2) of the 
2x2 block matrix with non- vanishing xi, i.e. those vec- 
tors which evolve perturbatively from those of A. There 
are additional eigenvectors of the form x = (0,^2), which 
can, however, be omitted, since their eigenvalues are of 
order 0{V) instead of 0{t). 

3.2. Solution of the non-linear eigenvalue problem 

In the numerical projection technique, to be outlined 
below, the second partition will consist of a single sector. 
Hence B will be the sub-matrix of the Hamiltonian cor- 
responding to basis vectors of that particular sector and 
it will consist of a kinetic term and an interaction term, 
B = B-\-By^ with By = Nf^VI^ where a I is a unit matrix. 
By virtue of this structure and the fact that \V\ the 
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numerical solution of the non-linear problem can be sim- 
plified significantly. We can expand the inverse in powers 
of By^ or rather {By — A)~^. The latter step is justified, 
since the lowest eigenvalues are of order t plus corrections 
of order hence By — X = 0{V). The expansion yields 



1 



1 



B-X 



Bv-X^B 



1 



Bv-X {By - A)2 
_ 1 



B 



E 



=0 



{Bv - A)-+i 



{-BY 



and the calculation of S\ becomes 



Sx = A + J2iX - Sv/)"^"+'^ EB'' . 



(4) 



u=0 



Note that the expressions (A — are only num- 
bers. For an iterative solution of the non-linear eigenvalue 
problem, has to be computed repeatedly for different 
values of A, but only the pre-factors (A — By) 



-(^+1) 



are 



modified. The matrices EB^ E\ however, are independent 
of A. Therefore, they can be calculated once and stored, 
since they are in general of moderate size, as will be dis- 
cussed later on. 

This expansion can also be used in Eq. |2] in order to 
evaluate the second part of the eigenvector, i.e. 



X2 



-{B - X)-^E^ = - ^(A - By)-^^+^^ B^ E^ 



The non-linear eigenvalue problem can now be solved 
iteratively, by using an initial approximation for A to build 
aSa, then solve the eigenvalue problem by suitable means 
and use the resulting eigenvalue as a new approximation 
for A. Let xx be the normalized eigenvector of Sx to the 
lowest eigenvalue f (A). In the next recursion A^"®^^ = 
£{X) is used as new value for the parameter A. 

A more sophisticated approach to obtain an improved 
value for A is provided by the Newton-Raphson method. 
We are actually seeking the roots of ^{X) = ^(A) — A. The 
Newton-Raphson approach yields 



A 



j'(A) _ 1 
*'(A) ~ ~£-'(A) - 1 



^(A) 



g'(A) 
^'(A) - 1 



A . 



The expression shows that this can be transformed to an 
equation of form A^"^^) = a£{\) + (1 - a)\ with 

1 



£'{\) - 1 

Note that in the simple iteration scheme a = 1. In order to 
calculate ^'(A) we can use the Hellmann-Feynman theorem 
resulting in 

£'{\) = x\S'xXx . 



The derivative of 5'a (Eq. [s]) with respect to A therefore 
reads 

51 = -E{B - A)-2^t ^ (5) 
which, like Eq. [4j can be expressed in a Taylor expansion 

CO 

S\X) = - ^{ly + 1)(A - By)-^^+^^x\EB^E^xx • (6) 

As before, the stored matrices EB^E^ can be reused. 

3.3. Iterative inclusion of higher sectors 

Now we are in the position to describe the numerical 
projection technique. To begin with, we consider the fol- 
lowing two partitions: the first partition may contain the 
lowest sectors = 0,1,..., iV^ The second partition 
consists of merely one sector, namely = -\- 1. The 
sizes of the partitions are pi and p2, respectively. The 
projection step, discussed before, is used to combine the 
two partitions. This step can now be invoked repeatedly 
in order to include increasingly higher sectors. 

a) We begin with the first partition, solve the eigenvalue 
problem for Ai and determine the corresponding uni- 
tary matrix Vi containing the corresponding pi eigen- 
vectors. 

b) The matrix for the first two partitions has the follow- 
ing block structure 



^2 = I 



Ai El 
El Bi 



c) Next we perform a unitary transformation on A2 pro- 
vided by the unitary matrix diag(Vi,/) resulting in 



^'-[eIv, b. 



(7) 



with Di := F^^AiVi being a pi-dimensional diagonal 
matrix. 

d) Now we solve the eigenvalue problem for 
A2V2 = V2D2 as outlined in the previous sec- 
tion and construct the eigenvector matrix for the 
entire matrix A2: 



Vo 



with / being a identity matrix of appropriate size. As 
discussed before V2 only contains the lowest pi eigen- 
vectors. Consequently, D2 is again a pi -dimensional 
diagonal matrix, 
e) Finally we include the blocks E2 and B2 of the next 
sector and perform again a unitary transformation 
with diag(y2, 1) as in Eq. [t] leading to 



A. 



D2 



V^E^ 
B2 



The recursion steps (d) and (e) are continued until the 
desired accuracy is reached. 
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Figure 2: Truncation procedure for the scheme presented in the text. 
Each partition is consecutively projected onto the original effective 
matrix and diagonalized. 

3.4' Truncation of the first partition 

The advantage of the numerical projection technique so 
far is the fact that (non-hnear) eigenvalue problems have to 
be solved for matrices of a size given by the first partition 
Pi, which is obviously smaller then the original size. 

It may happen that the first partition, even if it covers 
merely the first sector, is still very large. This is the case 
in the Hubbard model, where already the lowest sector for 
= yields pi = (^^)(^^^^), which at half-filling be- 
comes pi = (^/2)* ^^^^ we will try and reduce the size 
of the first partition. Firstly, we can exploit translational 
invariance and solve the eigenvalue problem in the individ- 
ual /c-space sectors. We can furthermore restrict the num- 
ber of vectors, which are retained in the first partition, 
based on some suitable criterion as to there importance 
for the lowest eigenvalue. We will use only the most ob- 
vious and simple criterion, namely to keep those vectors, 
which belong to the lowest eigenvalues D2 of Ai. More 
sophisticated criteria as in DMRG are also conceivable. 

We will use a predefined value pi < pi for the trunca- 
tion size of the first partition. That means, we solve the 
eigenvalue problem of Ai of size pi x pi and retain only 
those pi eigenvectors in Vi which correspond to the low- 
est eigenvalues. According to the recursion procedure, in 
the following recursion steps the diagonal matrices enter- 
ing the left-upper block Dn will all have the size pi and 
henceforth only (non-linear) eigenvalue problems of this 
size occur. 

4. Numerical experiments 

We begin with the spinless fermion model with periodic 
boundary condition (pbc) at half-filling. Due to the con- 
servation of the total momentum k only those basis vectors 
with the same k need to be taken into account. For the 
spinless fermion model at half-filling the ground state is 
degenerate and has k = We will therefore only take 
basis vectors into account, that have a total momentum 

= ± 1^ The selection of the total momentum leads to a 
significant reduction of the computational effort. On the 
one hand, the size pi of the initial partition is reduced 
roughly by L and on the other hand it can be avoided 



^We keep both, to allow real valued matrices and eigenvectors. 



that in the truncation step vectors are retained, that are 
irrelevant for the ground state of the original system. 

To begin with, we choose as first partition the first sec- 
tor with = 0, which consists of merely 2 states with 
particles on every other site (charge density wave), i.e. 
Pi = Pi =2. These two vectors are the real valued linear 
combinations of the basis vectors. 



L 


M 


^NPT [s] 


tAR [s] 






16 


12870 


0.0 


0.2 


-1.57 


0.007 


18 


48620 


0.0 


0.8 


-1.76 


0.010 


20 


184756 


0.1 


4.1 


-1.95 


0.012 


22 


705432 


0.3 


10.4 


-2.14 


0.014 


24 


2704156 


0.6 


87.1 


-2.33 


0.014 



Table 1: Run-time and accuracy comparison between NPT (^npt? 
^npt) and traditional Arnoldi eigensolver (ARPACK, ^ar, ^0 ) for 
the spinless fermion model = 10) at half-filling for different 

system sizes and pbc. The recursion was stopped at an absolute 
accuracy of 10~^. The matrices have size M x M. All computations 
in this paper have been performed on a computer with Intel Core 2 
Duo 3GHz, 2GB RAM. 

Tab. [1] exhibits a comparison of the run time for various 
system sizes. It can be seen that NPT is significantly faster 
than the highly optimized ARPACK routine (Arnoldi im- 
plementation) for increasing system size. We observe that 
even for 10^ x 10^-matrices the CPU time is less than 1 
sec. and more than a hundred times faster than the highly 
optimized ARPACK code. 

4.1. Dependence on the size of the initial partition 

In Fig. |3]the upper curve (green squares) corresponds to 
the case that the initial partition consists of the first sector 
{N^ = 0) only, with Pi = Pi = 2 and shows the depen- 
dence of the lowest eigenvalue on the maximum number of 
sectors, specified by N^^^^ included in the recursion. We 
observe that the recursion rapidly levels off at a moderate 
relative accuracy of 1 percent. The reason is due to the 
fact, that the vectors that are omitted during the recursion 
procedure, are relevant for a higher accuracy. 

One way to increase the accuracy is to choose a greater 
initial partition upon including all sectors up to iV^^S but 
still using a truncation size of the first partition pi = 2. 
The results are also depicted in Fig. |3] The turquoise 
upward triangles, e.g., depict the results for N^^^ = 2, i.e. 
the initial partition contains sectors up to = 2 and 
the subsequent recursions mix in A^^ = 3, . . . , N^^^. The 
relative accuracy increases by one order of magnitude but 
is still restricted. In all curves we observe that the result 
is converged when two more sectors are included beyond 
the ones present in the initial partition, which is obviously 
a consequence of the small truncation size pi = 2. The 
blue solid line indicates the result of the lowest eigenvalue 
obtained in the subspace spanned by the vectors of the 
sectors = 0, . . . , A^^^^. 
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Figure 3: Convergence of the lowest eigenvalue with the largest sector 
(A/"™^^) included in the recursion, depending on the initial partition. 
In all cases, pi = 2 was used, i.e. the 2 lowest eigenvectors of the 
initial partition where kept. The blue solid line indicates the exact 
values corresponding the various subspaces (see text). The results 
are compared with the exact eigenvalue Eq of the original matrix on 
a logarithmic scale. System parameters are L = 20, = 10, \ V/t\ = 
10. 



5. Three-partition Projection 

As discussed in the introduction, the central goal of the 
projection techniques is to describe the low- lying generic 
physics of the original quantum many-body system. We 
have seen, that NPT as opposed to standard projection 
techniques resulting in effective models, even yields quan- 
titatively good agreement with the result of the the orig- 
inal model. If still higher accuracy is the ultimate goal, 
then one can do even better then what has been discussed 
so far. It is possible to include in each recursion step two 
sectors at a time with little more CPU effort. To this end 
we consider a matrix which has 3x3 tridiagonal block 
structure, that results, if we add in each recursion (pro- 
jection) step two additional sectors to the first partition. 
The matrix has the structure 

A E 0\ 

B f \ . 
Ft C) 



4-2. Dependence on the truncation size 

So far we have seen that the accuracy is limited when 
the truncation size is restricted to 2, even if the size of the 
initial partition is increased and all sectors are iteratively 
included. One is therefore driven to study the dependence 
on the truncation size. To this end we consider one par- 
ticular system {L = 20, N = 10,\V/t\ = 10) and choose as 
initial partition the lowest two sectors with = and 
= 1. The size of the initial partition is pi = 182. In 
Tab. [2] the energies are listed for different truncation sizes 
pi. In all cases NPT is iterated until an absolute accuracy 
of = 10~^ is reached. We observe a clear improvement 
as compared to the case pi = 2, although the accuracy is 
certainly limited by the fact, that only the first two sectors 
are included in the initial partition. We find that the accu- 
racy achieved by using three sectors can never be reached. 



Pi 


^NPT 


E]siPT — Eo 
Eo 


^iters. [^] 


2 


-1.9565 


0.012 


0.66 


4 


-1.9573 


0.011 


1.34 


6 


-1.9594 


0.010 


2.25 


10 


-1.9629 


0.008 


3.76 


14 


-1.9713 


0.004 


6.79 


16 


-1.9749 


0.003 


9.21 


3 sectors 
exact 


-1.9793 
-1.9800 


0.001 





Table 2: Dependence of the lowest eigenvalue on the truncation size 
PI for a system with L = 20, = 10, \V/t\ = 10, ex = 10"^. Here, 
the initial partition consists of the first two sectors. In addition the 
exact lowest eigenvalue is given for a) the sub-matrix including the 
first three sectors and b) the original matrix. 




Figure 4: Comparison of the 2-partition and 3-partition approach 
for two different system sizes. The convergence of the relative uncer- 
tainty of the lowest eigenvalue is shown versus the maximum number 
of sectors included. 

Here, an equivalent relation to Eq. [3] can be obtained: 

(^A-E{B-XI- F{C - A/)-^Ft)"^ £;t^ xi = Xxi . 

The contribution of the matrix C shall be expanded up 
to the linear inverse term: 

A- xi-E{B - xiy^ E"^ 

-E{B - XI)-^F{C - XI)-^F\B - XI)-^E^ 

Again, the expression is expanded in terms of the kinetic 
part: 
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7. Conclusions 



E{B - \I)-^F{C - \I)-^F\B - Xiy^E"^ = 
1 1 1 



E 



A)^+i {By - A)/^+i {Cv - A)^+i 

e{-byf{-cyf\-bye^ . (8) 



Using this extension, the numerical result can be im- 
proved by orders of magnitude as is shown in Fig. [4j The 
results are for the spinless fermion model at strong cou- 
pling \V/t\ — 10 and half-filling. In both case the initial 
partition is given by the lowest sector without further trun- 
cation, i.e. V\—V\— 2. 

6. Results for the Hubbard model 

Next we present some results for the Hubbard model. 
Standard projection technique leads in the case of half- 





system 




model 




E-Eo 
Eo 




L = 
U = 


8, N^=N^ = 
10 


= 4, 


Hubbard 
NPT (pi = 20) 
Heisenberg 


-2.1767 
-2.1526 
-2.2604 


0.0110 
0.0385 


L = 
U 


10, = 
10 


= 5, 


Hubbard 
NPT (pi = 26) 
Heisenberg 


-2.7037 
-2.6572 
-2.8062 


0.0172 
0.0379 


1. = 

U 


10, = 
10 


= 4, 


Hubbard 
NPT (pi = 80) 
tJ 


-5.6698 
-5.7303 
-5.5282 


0.0107 
0.0250 


1. = 

U = 


10, = 
5 


= 5, 


Hubbard 
NPT (pi = 26) 


-4.9334 
-4.6032 


0.0669 






Heisenberg 


-5.6123 


0.1376 



Table 3: Comparison of the lowest eigenvalue obtained by NPT with 
those of the Heisenberg or t J-model and the exact result for the 
original Hubbard model. Periodic boundary conditions are assumed. 
Details of the NPT parameters are given in the text. One sectors 
were taken as initial partition (N^^^ = 1). 

filling to the Heisenberg model and away from half-filling 
to the t J- model [9 . These models correspond in NPT to 
the situation that the initial partition is the first sector 
with A^^ = (no double occupancies) and only the second 
sector is mixed in the projection step. On top of that, 
the non-linear eigenvalue problem is replaced by the linear 
one with A = in Eq. [2] In Tab. [3] the lowest eigenvalues 
for different system sizes and particle numbers are given. 
In the half-filled case the NPT results are compared with 
the exact results for the corresponding Heisenberg model 
and away from half-filling with those of the t J-model. In 
NPT we use the first sector {N^ = 0) as initial partition, 
restricted to /c = 0. For the truncation size in the projec- 
tion steps we used pi = 2. The NPT recursion steps are 
stopped at an absolute accuracy ex = 10~^. 

We compare NPT results and those of the effective mod- 
els with the exact lowest eigenvalue of the Hubbard model. 
We see that even in the strong coupling case {\U/t\ = 10), 
investigated here, NPT yields quantitatively better results 
than the effective models. 



The numerical projection technique (NPT), presented in 
this paper, is a numerical generalization of standard pro- 
jection technique, routinely used to derive effective mod- 
els for qualitative description of the generic low-energy 
physics of strongly correlated systems. The standard pro- 
jection technique corresponds to an approximate evalua- 
tion of second order perturbation theory in the inverse 
interaction parameter. 

NPT on the one hand a systematic application of the 
projection techniques to any many-body Hamiltonian, ir- 
respective of the complexity of the original model, also in 
cases where the analytic approach would lead to unman- 
ageable effective models. On the other hand it is also ap- 
plicable for moderate coupling parameters, where second 
order perturbation theory is not reliable enough. 

We have demonstrated that NPT yields fairly accurate 
results with little CPU-effort as compared to state-of-the- 
art exact diagonalization. Furthermore, by solving the 
non-linear eigenvalue problem and including higher order 
terms the approach yields significantly better results than 
traditional projection technique. 
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